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1  Grant  summary 

This  report  summarizes  the  research  completed  under  the  grant  “Design  of  Hierarchical  Hybrid 
Control  Systems”,  under  the  DARPA  Software  Enabled  Control  (SEC)  program,  administered  by 
AFRL  under  contract  F33615-99-C-3014,  which  was  awarded  to  Stanford  University  over  the  period 
10/15/99  to  01/31/05.  The  total  amount  awarded  was  $3,048,809.  The  Principle  Investigators  for 
this  grant  were  Professors  Claire  Tomlin  (PI)  and  Stephen  Boyd  (Co-PI). 

The  four  pillars  of  the  Stanford  SEC  program  were: 

1.  design  and  development  of  innovative  mathematical  methods  for  control  of  networks  of  com¬ 
plex  systems; 

2.  design  and  development  of  a  computational  tool  for  (close  to)  real-time  verification  of  complex, 
hybrid  systems; 

3.  design  and  development  of  the  embedded  avionics  for  each  aircraft  in  a  fleet  of  two  small  Un¬ 
manned  Air  Vehicles  (the  Stanford  DragonFly  UAVs),  and  the  coordination  and  cooperation 
of  this  development  with  the  Boeing  OCP  developers; 

4.  the  implementation  and  testing  of  the  hybrid  system  reachability  tool,  and  the  control  algo¬ 
rithms  for  networks  of  complex  systems,  on  the  Stanford  DragonFly  UAVs,  and  on  the  Boeing 
Capstone  Demonstration  Platform. 

2  Personnel 

Over  the  course  of  the  grant,  the  personnel  supported  were: 

‘Associate  Professor,  Aeronautics  and  Astronautics,  and,  by  courtesy,  of  Electrical  Engineering 
tomlin@stanford.edu,  (phone)  650-723-5164,  (fax)  650-723-3738,  Durand  250,  Stanford  CA,  94305-4035. 
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Faculty:  Claire  Tomlin  (Aero/Astro),  Stephen  Boyd  (EE),  David  Dill  (CS),  and  Zohar  Manna 
(CS). 

Postdoctoral  Research  Associates:  Haitham  Hindi  (EE,  now  at  PARC,  Palo  Alto);  Dusan 
Stipanovic  (Aero/ Astro,  now  Assistant  Professor,  University  of  Illinois  at  Urbana  Champaign); 
Henny  Siprna  (CS);  Mikael  Johansson  (EE,  now  Assistant  Professor,  KTH  Stockholm);  Jung  Soon 
Jang  (Aero/ Astro). 

Graduate  students:  Ian  Mitchell  (now  Assistant  Professor,  UBC);  Michael  Colon  (CS);  Gokhan 
Inalhan  (Aero/Astro,  now  postdoctoral  researcher,  MIT);  Inseok  Hwang  (Aero/Astro,  now  Assis¬ 
tant  Professor,  Purdue);  Jung  Soon  Jang  (Aero/Astro);  Meeko  Oishi  (ME,  now  postdoctoral  re¬ 
searcher,  NEON);  Lin  Xiao  (Aero/ Astro,  now  postdoctoral  researcher,  Caltech);  Alexandre  Bayen 
(Aero/Astro,  now  Assistant  Professor,  Berkeley). 


3  Research  Area  1:  Controlling  systems  over  networks 

The  continuing  drop  in  prices  for  data  networks  and  computers  has  created  a  strong  drive  to  let 
control  system  components  communicate  over  a  shared  communication  link  rather  than  using  sep¬ 
arate  dedicated  communication  channels.  However,  the  presence  of  a  communication  network  in 
the  control  loop  complicates  the  control  design.  Communication  over  a  shared  medium  introduces 
varying  time  delays,  occasional  loss  of  data  and  even  the  possibility  of  link  failures.  Traditional 
control  theory  gives  little  support  for  dealing  with  the  these  uncertainties,  and  an  understanding 
of  the  issues  in  networked  control  systems  has  just  recently  started  to  emerge.  So  far,  the  research 
on  networked  control  systems  has  been  exclusively  devoted  to  single-byte  data  transfer.  Today, 
however,  transfer  times  are  usually  a  negligible  part  of  the  communication  delay  and  most  network 
protocols  allow  us  to  send  several  bytes  of  data  at  no  extra  cost.  This  makes  it  interesting  to 
investigate  how  multiple-byte  transfers  can  be  exploited  to  improve  system  performance  and  ro¬ 
bustness  to  varying  time-delays,  occasional  loss  of  data,  and  possible  link  failures.  We  have  derived 
optimal  controllers  and  estimators  for  packet-oriented  communication,  and  shown  that  in  situations 
with  varying  time  delays,  they  perform  better  than  approaches  that  use  single-byte  data  transfer, 
without  incurring  any  additional  load  on  the  communication  network  [1], 

Now,  consider  jointly  allocating  communication  resources  and  designing  a  linear  control  system  in 
order  to  optimize  system  performance.  This  problem  is  in  general  not  convex,  but  can  be  solved 
heuristically  in  a  way  that  exploits  the  special  structure  of  the  communication  resource  allocation 
problems.  Numerical  examples  show  that  our  approach  of  jointly  optimizing  both  the  system  and 
communication  resources  can  have  very  significant  benefits  over  the  traditional  approach  [2], 

Finally,  having  established  the  SRRA  framework  with  theoretical  channel  capacity  functions  in 
previous  work,  we  have  moved  one  step  further  to  solve  the  SRRA  problem  with  practical  fading 
wireless  channels.  We  showed  that  the  joint  optimization  of  data  routing,  resource  allocation  and 
power  control  over  fading  states  can  be  efficiently  solved  by  convex  optimization  techniques,  which 
opens  the  door  for  practical  application  of  the  SRRA  framework.  We  also  extended  the  SRRA 
framework  to  heterogeneous  networks,  where  wired  and  wireless  networks  with  different  media 
access  schemes  are  connected  together  to  facilitate  ubiquitous  data  communication.  Our  SRRA 
framework  can  be  extended  naturally  to  handle  such  complex  systems. 
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4  Research  Area  2:  Computing  Reachable  Sets  for  Hybrid  Sys¬ 
tems 

The  hybrid  systems  framework  provides  an  appealing  means  for  verifying  the  safety  of  dynamical 
systems.  We  address  safety  verification  as  a  reachability  problem:  Given  an  unsafe  subset  of  the 
system  state  space  and  an  initial  state,  is  the  former  reachable  from  the  latter?  If  it  is  reachable 
despite  any  controllable  actions  that  we  can  take  then  the  system  is  provably  unsafe. 

The  continuous-time  nonlinear  dynamics  of  the  system  need  to  be  considered  in  assessing  reach¬ 
ability.  Our  formulation  requires  the  solution  of  a  Hamilton- Jacobi  partial  differential  equation, 
and  a  grid-based  numerical  solution  approach  based  on  level  set  methods  is  used  for  this  purpose. 
Simulation  examples  are  presented  for  three  flight  management  applications:  two-aircraft  colli¬ 
sion  avoidance,  the  related  problem  of  conflict  resolution,  and  ensuring  safety  during  final  landing 
approach. 

The  exact  reachability  computation  is  prey  to  the  curse  of  dimensionality:  its  computational  com¬ 
plexity  is  exponential  with  respect  to  the  continuous  dimension.  We  also  developed  alternative 
approaches,  which  are  based  on  over-approximating  the  reachable  set  of  states  with  a  polyhedron. 
This  is  also  computationally  intractable  since  the  propagation  of  the  system’s  dynamics  will  result 
in  a  potentially  unlimited  number  of  constraints  (faces  of  the  polyhedron),  but  we  have  developed 
a  novel  technique  for  identifying  and  pruning  redundant  and  irrelevant  constraints.  This  technique 
promises  to  be  computationally  feasible  for  very  high  dimensional  problems.  The  basis  of  the 
approach  is  the  computation  of  the  maximum  volume  ellipsoid  contained  in  a  polyhedron,  a  com¬ 
putation  that  can  be  formulated  as  a  convex  optimization  problem  (for  which  global  and  efficient 
algorithms  are  available). 

4.1  Introduction 

For  about  the  past  ten  years,  researchers  in  the  traditionally  distinct  fields  of  control  theory  and 
computer  science  verification  have  proposed  models,  and  verification  and  controller  synthesis  tech¬ 
niques  for  complex,  safety  critical  systems.  The  area  of  hybrid  systems  theory  studies  systems  which 
involve  the  interaction  of  discrete  event  and  continuous  time  dynamics,  and  provides  an  effective 
modeling  framework  for  complex  continuous  systems  with  large  numbers  of  operating  modes.  Ex¬ 
amples  include  continuous  systems  controlled  by  a  discrete  logic  such  as  the  autopilot  modes  for 
controlling  an  aircraft,  systems  of  many  interacting  processes  such  as  air  or  ground  transportation 
systems  in  which  discrete  dynamics  are  used  to  model  the  coordination  protocols  among  processes, 
or  continuous  systems  which  have  a  phased  operation,  such  as  biological  cell  growth  and  division. 

The  current  and  potential  impact  of  hybrid  systems  lies  in  the  confluence  of  computational  methods 
from  control  theory  and  from  formal  methods  in  computer  science  verification.  In  the  examples 
mentioned  above,  the  system  dynamics  are  complex  enough  that  traditional  analysis  and  control 
methods  based  solely  on  differential  equations  are  not  computationally  feasible;  analysis  based 
solely  on  discrete  event  dynamics  ignores  critical  system  behavior.  Our  interest  lies  in  developing 
computational  tools  for  analyzing  and  controlling  the  behavior  of  hybrid  systems;  our  eventual 
goal  is  to  develop  a  real-time  tool  to  provide  online  verification  that  a  hybrid  system  satisfies  its 
specified  behavior.  This  goal  addresses  one  of  the  central  themes  of  Software  Enabled  Control,  in 
which  we  seek  a  systematic  methodology  for  the  design,  validation,  and  implementation  of  control 
software. 


3 


In  our  work,  the  system  specification  that  we  are  most  interested  in  is  that  of  system  safety,  which 
asks  the  question:  Is  a  potentially  unsafe  configuration  of  the  system  reachable  from  an  initial 
configuration?  More  important  for  control  theory,  given  a  set  of  desired  configurations,  it  is  crucial 
to  be  able  to  design  the  hybrid  system  control  inputs  to  achieve  these  configurations.  Previous  work 
in  this  area  had  focused  on  hybrid  systems  with  very  simple  continuous  dynamic  equations  (such  as 
clocks,  or  linear  decoupled  systems).  Our  research  has  three  components:  the  problem  formulation 
for  computing  exact  reachable  sets  of  hybrid  systems  and  the  design  of  a  software  tool  to  perform 
such  calculations;  the  development  of  an  algorithm  for  computing  over  approximations  of  reachable 
sets  which  works  efficiently  in  high  dimensions;  the  design  of  a  real-time  aircraft  testbed  for  these 
algorithms.  We  will  focus  on  the  first  component  in  this  report,  we  briefly  outline  our  work  in  the 
second. 

4.2  Hybrid  System  Model 

The  main  difference  between  our  hybrid  system  model  and  leading  models  in  the  literature  (see  for 
example  the  work  of  Alur  and  Henzinger  [3]  on  linear  hybrid  automata)  is  that  we  include  accurate, 
nonlinear  models  of  the  hybrid  system  continuous  dynamics.  The  model  and  algorithm  have  been 
developed  jointly  with  Lygeros  and  Sastry,  full  details  are  presented  in  [4], 

A  hybrid  automaton  is  a  finite  state  machine  with  discrete  states  {91,92,  ,Qk},  in  which  each 

discrete  state  has  associated  continuous  dynamics  x  =  f(qi,x,v),  with  x  £  Rn,  and  continuous 
inputs  and  disturbances  v  =  ( u,d ),  where  u  £  1/  and  d  £  D. 

Definition  1  (Hybrid  Automaton).  A  hybrid  automaton  H  is  a  collection 

H  =  (■ Q ,  X,  £,  V,  Init,  /,  Inv,  R) 


where 


•  Q  Li  X  is  the  set  of  state  variables,  with  Q  a  finite  set  of  discrete  states,  and  X  =  Mn; 

•  S  =  Si  U  £2  is  a  finite  collection  of  discrete  input  variables,  where  Si  is  the  set  of  discrete 
control  inputs,  and  S2  is  the  set  of  discrete  disturbance  inputs; 

•  V  =  U  U  D  is  the  set  of  continuous  input  variables,  where  U  is  the  set  of  continuous  control 
inputs,  and  D  is  the  set  of  continuous  disturbance  or  uncontrollable  inputs;  we  denote  the 
spaces  of  input  (disturbance)  trajectories  as  the  sets  of  piecewise  continuous  functions  U  ( V 
respectively) ,  which  take  values  in  U  (D  respectively) ; 

•  Init  C  Q  x  X  is  a  set  of  initial  states; 

•  f:  QxXxV^TX  is  a  vector  field  describing  the  evolution  of  x  for  each  q  £  Q;  f  is 
assumed  to  be  globally  Lipschitz  in  X  (for  fixed  q  £  Q)  and  continuous  in  V; 

•  Inv  CQxXxHxV  is  called  an  invariant,  and  defines  combinations  of  states  and  inputs 
for  which  continuous  evolution  is  allowed; 

•  R:Qx  X  x  SxF->  2®xX  is  a  reset  relation,  which  encodes  the  discrete  transitions  of  the 
hybrid  automaton. 
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We  refer  to  (q,  x)  E  Q  x  X  as  the  state  of  H  and  to  ((<ti,  (T2),  (■ u ,  d))  Sflxhas  the  input  of  H.  The 
control  actions  (a\,u)  model  those  inputs  over  which  the  designer  has  control,  and  the  disturbance 
actions  (<72,  d)  model  inputs  over  which  the  designer  has  no  control,  such  as  uncertainties  in  the 
actions  or  behaviors  of  the  other  aircraft  in  the  system.  We  assume  that  the  designer  has  complete 
knowledge  of  the  bounds  on  these  disturbance  actions  (£2 , -D). 

Associated  to  the  hybrid  automaton  H  is  a  specification  which  describes  the  condition  one  would 
like  the  system  to  satisfy.  Our  work  has  been  motivated  by  verification  and  synthesis  for  safety 
critical  applications,  and  as  such  we  have  been  primarily  interested  in  safety  specifications.  These 
specifications  are  encoded  as  subsets  of  the  state  space  of  the  hybrid  system:  the  set  Go  Cl  Q  x  X 
is  that  subset  which  is  defined  a  priori  to  be  unsafe.  The  problem  specification  is  therefore  to  keep 
the  system  state  outside  of  the  unsafe  set. 

4.3  Exact  Reach  Set  Computation  using  Level  Sets 

In  this  section,  we  describe  the  development  of  a  general  purpose  tool  for  exact  reachable  set 
computation — the  core  of  which  is  a  new  variant  of  a  “local  level  set”  algorithm  that  efficiently 
computes  an  accurate  representation  of  the  reachable  set  boundary.  We  demonstrate  on  an  example 
the  numerical  convergence  of  our  computation  by  analyzing  the  results  as  the  continuous  state  space 
grid  is  made  finer,  a  standard  method  of  validation  for  scientific  computing  codes.  In  this  way,  we 
show  that  high  accuracy  can  be  achieved  at  the  cost  of  increased  computational  time  and  space. 
We  illustrate  our  tool  on  a  single  mode  aircraft  conflict  resolution  example  [5,  6],  a  three  mode 
conflict  resolution  example,  as  well  as  an  example  of  a  six  mode  commercial  aircraft  auto-lander 
[7],  which  exhibits  nondeterminism  and  cycles  in  its  discrete  behavior. 

Our  motivation  for  this  component  of  the  research  stems  from  the  belief  that  for  many  applications 
of  hybrid  systems,  it  is  important  to  be  able  to  accurately  represent  the  reachable  set.  We  have 
dealt  primarily  in  the  safety  verification  of  avionic  systems,  where  accurate  representation  of  the 
safe  region  of  operation  translates  into  the  ability  to  operate  the  system  closer  to  the  boundaries 
of  that  region,  at  a  higher  performance  level  than  previously  allowed.  For  very  high  dimensional 
state  spaces,  additional  logic  (such  as  projection  operators)  or  new  techniques  (such  as  the  convex 
overapproximations  of  the  next  section)  will  be  needed;  however,  our  results  in  this  paper  show  that 
it  is  feasible  to  do  exacting  computation  for  hybrid  systems  with  nonlinear  continuous  dynamics 
in  three  continuous  state  dimensions  and  six  discrete  modes,  and  we  believe  it  will  be  feasible  to 
extend  this  up  to  five  continuous  dimensions  and  large  numbers  of  discrete  modes. 

4.4  Reachability  for  Hybrid  Systems 

In  this  section  we  summarize  the  general  framework  for  handling  the  interaction  between  discrete 
and  continuous  dynamics  (following  [4]). 

Fundamentally,  reachability  analysis  in  discrete,  continuous  or  hybrid  systems  seeks  to  partition 
states  into  two  categories:  those  that  can  reach  a  given  target  set,  and  those  that  cannot.  We  will 
label  these  two  sets  of  states  G  and  E  =  Gc  respectively  (the  symbol  G  has  been  used  traditionally  to 
represent  unsafe  sets,  E  is  used  as  a  symbol  for  “escape”  sets).  Any  inputs  to  the  hybrid  automata 
are  assumed  to  lie  in  bounded  sets  and  to  have  the  goal  of  locally  maximizing  or  minimizing 
the  reachable  set:  at  each  iteration,  the  reachability  algorithm  chooses  values  for  inputs  that 
maximize  the  size  of  G  and  values  for  inputs  £e  that  minimize  the  size  of  G  (and  hence  maximize  the 
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mode 
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•  •  • 


•  •  • 


•  •  • 

•  •  • 


K 


iteration 

0 


(mode,  iteration)  =  (j,  i) 


Figure  1:  Iterative  Reachability  Algorithm:  Showing  detail  of  iteration  for  discrete  mode  j  at 
iteration  i. 


size  of  E).  Any  nondeterminism  in  the  transition  relation  is  also  utilized  to  consistently  maximize 
or  minimize  G,  depending  on  the  goal  of  the  reachability  computation.  For  hybrid  automata,  the 
discrete  inputs  a  and  continuous  inputs  v  can  be  assigned  to  the  two  categories  £g  =  (ctg,^g)  and 

=  (o'EjJ'e)  according  to  whether  they  seek  to  maximize  or  minimize  G.  For  example,  for  our 
hybrid  system  model  with  control  actions  (ay  ,u)  and  disturbance  actions  (oy  ,d),  we  would  have 
that  £g  =  (og^g)  =  {<t2,  d),  and  £E  =  (oe^e)  =  {<7\,u). 

The  reachability  computation  follows  an  iterative,  two  stage  algorithm  shown  graphically  in  Fig¬ 
ure  1.  The  outer  iteration  computes  reachability  over  the  discrete  switches,  producing  iterates  G^ 
and  Ej  at  iteration  i  =  0,-1,  —2, . . . ,  where  a  negative  index  is  used  to  indicate  that  the  algorithm 
is  initialized  with  a  target  set  and  run  backward  to  compute  those  states  which  can  reach  the  target 
set.  The  inner  iteration  runs  a  separate  continuous  reachability  problem  in  each  of  the  discrete 
modes  j  =  1,  2, ...  K  to  compute  the  estimates  Gj  and  Ej.  We  define  the  “switch”  sets 

•  Gj  contains  all  states  in  mode  j  from  which  a  discrete  transition  to  a  state  in  Gj_i  (typically 
a  state  in  another  mode)  can  be  forced  to  occur  through  the  application  of  a  discrete  input 
<jq\  these  states  will  be  defined  by  the  invariant  of  mode  j  and  the  guards  of  the  transitions 
from  mode  j. 

•  Ej  contains  all  states  from  which  a  discrete  transition  to  a  state  in  E?;_i  can  be  forced  to  occur 
through  the  application  of  a  discrete  input  uE;  these  states  are  also  defined  by  the  invariant 
of  mode  j  and  the  guards  of  transitions  from  mode  j. 

Then  the  goal  of  the  continuous  reachability  tool  is  to  identify  the  “flow”  sets 

•  Gj  ( t )  contains  states  from  which  for  all  z^E  there  exists  vq  that  will  force  the  resulting  trajectory 
to  flow  into  Gj_,  U  Gj  within  time  t. 

•  E j{t)  contains  states  from  which  there  exists  z^E  that  for  all  uq  will  force  the  resulting  trajec¬ 
tories  to  flow  into  E j  within  time  t  or  to  stay  outside  of  G-(_1  U  Gj  for  at  least  time  t. 

Note  that  in  some  problems  the  order  of  the  existential  and  universal  quantifiers  in  the  definition 
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above  must  be  reversed.  Given  these  sets 


K 

G i  =  M  Gji  G  =  lim  Gj, 

^  i— >oo 

3= 1 

K 

E?:  =  M  Ej,  E  =  lim  Eu 

4— XX) 

3= 1 

where  Gg  is  the  set  of  initial  conditions  of  the  reachability  problem  and  Eg  =  (Gg)c.  Simple 
modifications  of  this  algorithm  suffice  to  solve  finite  time  reachability  problems. 

The  procedure  described  above,  developed  in  [4,  6],  was  motivated  by  the  work  of  [8,  9]  for  reach¬ 
ability  computation  and  controller  synthesis  on  timed  automata,  and  that  of  [10]  for  controller 
synthesis  on  linear  hybrid  automata.  In  that  development  the  reachability  problem’s  objective  was 
to  determine  E  -the  largest  controllable  invariant  subset  of  the  state  space — by  computing  the  set 
of  states  G  which  were  reachable  in  backwards  time  from  the  set  of  predefined  unsafe  states.  In 
terms  of  the  definitions  above,  control  inputs  from  this  problem  lie  in  £e  and  disturbance  inputs  in 
£g-  For  safety,  any  model  nondeterminism  would  be  used  to  maximize  the  unsafe  set  G. 


G  ■  =  km  Gj(t), 

t — XX) 

Ej  =  lim  Ej(t), 

t — XX) 


4.5  Continuous  Reachability  using  Level  Sets 


While  practical  algorithms  for  computing  discrete  reachability  over  many  thousands  of  states  have 
been  designed  and  implemented,  determination  of  continuous  reachability  for  even  low  dimensional 
systems  is  still  an  open  problem.  The  continuous  portion  of  a  hybrid  reachability  problem  requires 
methods  of  performing  four  key  operations  on  sets:  unions,  intersections,  tests  of  equality,  and 
evolution  according  to  the  discrete  mode’s  continuous  flow  field.  The  choice  of  representation  for 
sets  dictates  the  complexity  and  accuracy  of  these  operations;  consequently,  continuous  reachability 
algorithms  can  be  classified  according  to  how  they  represent  sets. 

For  our  exact  representation  scheme,  we  characterize  the  set  being  tracked  implicitly  by  defining 
a  “level  set  function”  J(x,t)  throughout  the  continuous  state  space  which  is  negative  inside  the 
set,  zero  on  its  boundary,  and  positive  outside,  and  which  encodes  the  initial  data  in  J(x,0).  The 
intersection  of  two  such  sets  is  simply  the  maximum  of  their  level  set  functions  at  each  point  in  state 
space,  and  the  union  is  the  minimum;  a  variety  of  easily  implemented  equality  tests  are  possible. 
Evolution  of  a  level  set  under  a  nonlinear  flow  field  is  governed  by  the  Hamilton-Jacobi  (HJ)  partial 
differential  equation  (PDE)  (see,  for  example,  [5]) 


8J(x,  t) 
dt 


min{0,  max  min  f(x,  vm  in,  J(x,  f)}, 

^min  ^max 

min{0,  H(x,  VJ(x,  t))} 


(1) 


where  zzmin  are  those  continuous  inputs  trying  to  minimize  the  size  of  the  set  being  tracked,  and 
rmax  are  those  inputs  trying  to  maximize  its  size.  The  order  of  the  optimization  must  be  chosen 
appropriately  for  the  situation.  The  first  minimization  on  the  right  hand  side  of  (1)  ensures  that 
the  set  can  only  grow  (states  that  are  unsafe  cannot  become  safe  at  a  future  time).  That  the 
subzero  level  sets  of  the  viscosity  solution  of  (1)  yields  the  correct  reachable  set,  is  proved  in 
[11].  The  implicit  representation  has  a  number  of  advantages  when  compared  with  the  explicit 
representations  that  other  researchers  are  pursuing,  including  a  conceptually  simple  representation 
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of  very  general  sets  and  a  size  which  is  independent  of  the  complexity  of  the  set  (although  it 
grows  exponentially  with  dimension).  In  addition,  a  set  of  sophisticated  numerical  techniques  to 
accurately  solve  PDEs  may  be  drawn  upon  for  computation.  In  the  remainder  of  this  section, 
we  focus  on  the  representation  (1),  and  assume  that  the  modeler  can  compute  the  appropriate 
optimization  over  inputs  in  (1)  if  given  x  and  VJ(x,f). 

The  HJ  PDE  (1)  is  well  known  to  have  complex  behavior.  Even  with  smooth  initial  data  J(x,0) 
and  continuous  Hamiltonian  JJ(x,VJ),  the  solution  J(x,t )  can  develop  discontinuous  derivatives 
in  finite  time;  consequently,  classical  infinite  time  solutions  to  the  PDE  are  generally  not  possible. 
In  the  quest  for  a  unique  weak  solution  Crandall  and  Lions  introduced  the  concept  of  the  viscosity 
solution  [12].  For  most  problems  of  interest,  finding  the  analytic  viscosity  solution  is  not  possible 
(see  [13]  for  cases  in  which  it  is  possible),  and  so  we  seek  a  numerical  solution.  We  approximate  the 
solution  of  (1)  on  a  Cartesian  grid  of  nodes.  Three  terms  in  the  equation  must  be  approximated  at 
each  node,  based  on  the  values  of  the  level  set  function  at  that  node  and  its  neighbors:  the  gradient 
VJ,  the  Hamiltonian  H,  and  the  time  derivative  .  We  discuss  each  of  these  separately. 

In  each  dimension  at  each  grid  point  there  exist  both  left  and  right  approximations  of  the  gradient 
VJ,  depending  on  which  neighboring  grid  points’  values  are  used  in  the  finite  difference  calculation. 
We  label  the  vector  of  left  approximations  VJ-,  the  vector  of  right  approximations  VJ+,  and  will 
see  below  that  VJ-,  VJ+  or  some  combination  of  the  two  will  be  used  to  compute  the  numerical 
Hamiltonian  H.  The  accuracy  of  a  derivative  approximation  is  measured  in  terms  of  the  order  of 
its  local  truncation  error;  an  order  p  method  has  error  1 1 V J  —  VJ±||  =  0( Axp).  At  the  current 
time,  we  have  implemented  the  basic  first  order  accurate  approximation  for  speed  and  a  weighted, 
essentially  non-oscillatory  fifth  order  accurate  approximation  for  high  fidelity.  We  have  chosen  to 
use  the  well  studied  Lax-Friedrichs  numerical  Hamiltonian  approximation  H  [14] 

H(x,  VJ-,  VJ+)  =  H(x,  VJ~+VJ+)  -  \aT{VJ+  -  VJ-),  (2) 

where  H{x,  VJ)  is  given  by  (1)  and  the  term  containing  the  vector  coefficient  a  is  a  high  order 
numerical  dissipation  added  to  damp  out  spurious  oscillations  in  the  solution.  The  time  derivative 
of  the  PDE  is  handled  by  the  method  of  lines:  the  value  of  the  level  set  function  J  at  each  node 
is  treated  as  an  ODE  ^  =  H,  with  H  given  by  (2).  General  ODE  solvers,  such  as  Runge-Kutta 
(RK)  schemes,  can  then  be  applied. 

The  Hamilton- Jacobi  equation  (1)  describes  the  evolution  of  the  level  set  function  over  all  of  space. 
But  we  are  only  interested  in  its  zero  level  set;  thus,  we  can  restrict  our  computational  updates 
to  nodes  near  the  boundary  between  positive  and  negative  J(x,t ) — an  idea  variously  called  “local 
level  sets”  [15]  or  “narrowbanding”  [16].  We  have  implemented  a  new  variant  of  this  method  in  our 
code.  Because  the  boundary  is  of  one  dimension  less  than  the  state  space,  considerable  savings  are 
available  for  two  and  three  dimensional  problems.  If  the  number  of  nodes  in  each  dimension  is  n 
(proportional  to  Ax-1)  and  the  dimension  J,  the  total  number  of  nodes  is  0(nd)\  with  the  restriction 
on  the  timestep  necessary  for  numerical  stability,  the  total  computational  cost  is  0(nd+1).  With 
local  level  sets,  we  reduce  computational  costs  back  down  to  0(nd). 

4.6  Examples 

In  this  section  we  present  three  examples:  a  validation  of  our  numerical  implementation  on  a 
single  mode,  three  dimensional  aircraft  collision  avoidance  example  (see  [6,  5]  for  details),  and 
two  multimodal  examples.  Once  a  method  of  determining  continuous  reachability  is  available,  the 


discrete  iteration  of  the  algorithm  described  in  Section  4.4  is  relatively  straightforward.  In  fact,  for 
discrete  transition  graphs  with  no  cycles  it  is  possible  to  order  the  continuous  reachability  problems 
such  that  no  discrete  iteration  is  required  (the  three  mode  example  below) .  In  order  to  examine  the 
complications  induced  by  discrete  cycles — such  as  how  to  avoid  zenoness,  in  what  order  to  execute 
the  continuous  reachability  problems,  and  how  to  determine  which  switches  are  active —  we  present 
an  example  representing  the  landing  of  a  civilian  airliner. 

Numerical  Validation  of  Aircraft  Collision  Avoidance 

This  example  features  a  control  aircraft  trying  to  avoid  collision  with  a  disturbance  aircraft,  where 
both  aircraft  have  fixed  and  equal  altitude,  speed  and  turning  radius — they  may  only  choose  which 
direction  they  will  turn: 

xr  =  —vu  +  Vd  cos  tp r  +  uyr ,  yr  =  Vd  sin  ipr  —  uxr,  ipr  =  d  —  u, 

where  vu  =  Vd  =  5  are  the  aircraft  speeds,  xr  and  yr  are  the  relative  planar  location  of  the  aircraft 
and  ipr  is  their  relative  heading.  The  inputs  |u|  <  1  and  |d|  <  1  are  the  control’s  and  disturbance’s 
respective  turn  rates.  The  initial  unsafe  set  J(x,  0)  is  the  interior  of  the  radius  five  cylinder  centered 
on  the  tpr  axis.  Choosing  optimal  inputs  according  to  (1)  with  vq  =  =  d  and  zzE  =  ;ymin  =  u, 

we  get  the  optimal  Hamiltonian: 

H(x,p )  =  -piVu  +  Pivd  COS  Ipr  +  p2Vd  silllpr  +  \piyr  -p2Xr  P3 1  \P3  \  ■ 

Using  our  C++  implementation,  grid  sizes  corresponding  to  50,  70,  100,  140,  and  200  nodes  in 
each  dimension  were  tried  with  a  low  order  accurate  scheme  (first  order  space  and  time,  hereafter 
referred  to  as  the  “(1,1)”  scheme)  and  a  high  resolution  scheme  (fifth  order  space  and  second  order 
time,  hereafter  the  “(5,2)”  scheme).  On  the  eight  million  node  finest  grid — only  around  10%  of 
which  is  being  actively  updated  on  any  one  timestep  by  the  local  level  set  algorithm — execution 
time  for  the  (5,2)  scheme  was  about  eighteen  hours  on  a  Sun  UltraSparc  II  with  lots  of  memory. 
Reducing  the  grid  size  in  half  results  in  the  expected  eightfold  savings  in  memory  and  time;  hence, 
the  coarsest  grid  takes  only  fifteen  minutes  with  the  (5,2)  scheme. 

Results  are  visualized1  by  the  zero  level  isosurface  of  the  unsafe  reachable  set  G,  shown  in  Figure  2. 
On  the  left  is  a  head-on  view  of  the  (5,2)  solution.  On  the  right  is  a  zoomed  overhead  view  of  the 
point  of  the  bulge  computed  by  the  (1,1)  scheme  for  several  grid  sizes. 

We  compare  solutions  on  the  four  coarser  grids  to  the  solution  on  the  finest  grid,  using  linear 
interpolation  on  the  finest  grid  if  necessary.  Figure  3  demonstrates  that  the  scheme  is  converging 
to  the  finest  grid’s  solution  of  (1)  at  approximately  a  linear  rate  in  both  average  error  and  pointwise 
maximum  error.  We  cannot  expect  to  show  a  higher  order  convergence  rate  because  of  the  linear 
interpolation  used  to  evaluate  the  error. 

Two  conclusions  can  be  drawn  from  Figures  2  and  3.  First,  for  this  example,  low  order  schemes  are 
not  at  all  competitive  in  terms  of  accuracy  with  the  (5,2)  scheme:  our  (5,2)  implementation  can 
produce  more  accurate  results  in  about  fifteen  minutes  using  only  the  coarsest  grid.  Second,  the 
pointwise  maximum  error  of  the  (5,2)  scheme  is  always  less  than  the  grid  spacing,  so  if  a  5CU1  =  2% 
error  is  tolerable  for  this  application,  only  this  fastest,  coarsest  grid  need  ever  be  run. 

1  Figure  2  and  Figure  8  visualize  some  level  set  surfaces  as  triangular  meshes;  these  are  not  the  meshes  on  which 
the  Flamilton-Jacobi  PDE  was  solved,  but  rather  an  artifact  of  three  dimensional  Matlab  visualization  techniques. 


9 


(5,2)  scheme  on  1001  grid 


(1,1)  scheme  on  50d  (inner),  100J  (middle) 
and  2003  (outer)  grids 


>r  o 


Figure  2:  Reachable  Set  for  Aircraft  Collision  Avoidance  Example 


Three  Mode  Conflict  Resolution 

Consider  the  following  conflict  resolution  problem  between  two  aircraft  (i  E  {1,2}),  which  is  pre¬ 
sented  in  [17].  Let  (xr,yr,ipr)  E  M2  x  [— 7r,7r)  represent  the  relative  position  and  orientation  of 
aircraft  2  with  respect  to  aircraft  1.  In  terms  of  the  absolute  positions  and  orientations  of  the  two 
aircraft,  (xi,  yi,  ipi)  for  z  =  1,2,  it  may  be  verified  that  xr  =  cosip\{x2  —  x\)  +  sin '0i(j/2  —  Di),  Ur  = 
-sin'0i(ai2  -  x\)  +  cosi/’i(y2  -  yi),  ipr  =  i>2  -  Vh,  and  thus 

Xr  =  —Vi  +  V2  COS  1pr  +  UiUr 

yr  =  V2  sin  ipr  —  ui\xr  (3) 

tpr  =  <^2  -  Ob 


Figure  3:  Convergence  of  (5,2)  Scheme  to  Finest  Grid  Solution  (Jn  is  the  solution  J(x,t )  on  a  grid 
size  of  n) 
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Figure  4:  Two  aircraft  in  two  modes  of  operation:  in  the  first  and  third  figures,  the  aircraft  follow 
a  constant  heading  and  velocity  (in  Mode  1 );  in  the  second  figure,  the  aircraft  follow  a  half  circle 
(in  Mode  2). 


xr  =  -u  +dcosWr 
yr=  d  sin  Vr 
yr=0 
z  =  0 


xr  =  - u  +  d  cosxVr  +  yr 
yr=  d  sin^r-  xr 

%=0 
z  =1 


Xr  =  -u  +dcOSx\lr 
yr=  d  sin  Vr 
¥r=0 
z  =0 


Figure  5:  In  q\  both  aircraft  follow  a  straight  course,  in  qo  a  half  circle,  and  in  qs  both  aircraft 
return  to  a  straight  course. 


where  v%  is  the  velocity  along  the  body  axis  of  aircraft  i ,  and  u>i  is  the  rate  of  change  of  heading. 
A  conflict  occurs  when  a  pair  of  aircraft  incur  a  lateral  spacing  of  less  than  5  nautical  miles;  for 
safety  of  the  maneuver,  the  relative  position  ( xr,yr )  must  remain  outside  of  the  interior  of  the  5 
mile  radius  disk  centered  at  the  origin  {{xr,yr,i^r)  :  x2  +  y2  <  52}. 

Consider  the  maneuver  illustrated  in  Figure  4,  with  the  protocol  of  the  maneuver  defined  as  follows: 
the  aircraft  are  nominally  in  Mode  1,  in  which  the  aircraft  fly  straight  and  level  with  given  constant 
velocities  Vi  and  constant  heading  at  a  certain  relative  separation  distance  each  aircraft  turns  to 
its  right,  follows  an  arc  of  a  circle  in  Mode  2,  in  which  the  velocity  is  held  at  a  prescribed  constant 
value  Vi ,  until  it  intersects  its  original  trajectory,  then  turns  to  its  right  and  returns  to  its  desired 
trajectory.  The  model  allows  instantaneous  changes  in  heading  which  is  an  obvious  abstraction  of 
the  true  aircraft  dynamics  and  is  treated  here  in  order  to  clearly  illustrate  the  algorithm;  a  model 
with  heading  capture  dynamics  is  presented  in  [17].  In  each  mode,  the  continuous  dynamics  may 
be  expressed  in  terms  of  the  relative  motion  of  the  two  aircraft  (3):  in  Mode  1,  uii  =  0  for  i  =  1,2 
and  in  Mode  2,  cuj  =  1  for  i  =  1,2.  We  assume  that  both  aircraft  switch  modes  simultaneously,  so 
that  the  relative  orientation  ipr  is  constant.  This  assumption  simply  allows  us  to  display  the  state 
space  in  two  dimensions,  making  the  results  easier  to  present.  The  problem  statement  is  therefore: 
generate  the  relative  distance  between  aircraft  at  which  the  aircraft  may  switch  safely  from  Mode 
1  to  Mode  2,  and  a  turning  radius  r  in  Mode  2,  to  ensure  that  the  5  nautical  mile  separation  is 
maintained. 

The  dynamics  of  the  maneuver  can  be  encoded  by  the  hybrid  automaton  of  Figure  5,  where  q\ 
corresponds  to  Mode  1  before  the  maneuver,  q2  corresponds  to  Mode  2  “avoid  mode”,  and  q 3 
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corresponds  to  Mode  1  after  the  avoid  maneuver  has  been  completed.  There  is  one  discrete  control 
input  <7i,  such  that  the  switch  from  a\  =  0  to  =  1  triggers  the  transition  from  q\  to  qi .  The 
transition  from  q2  to  q3  is  required  to  take  place  after  the  aircraft  have  completed  a  half  circle: 
note  that  with  ujt  =  1,  for  i  =  1,2,  it  takes  n  time  units  to  complete  a  half  circle.  The  continuous 
state  space  is  augmented  with  a  timer  z  6  M+  to  force  this  transition.  Let  x  =  (xr,yr,ijjr,  z)T . 
At  each  transition,  both  aircraft  change  heading  instantaneously  by  it/ 2  radians;  we  represent  this 
with  the  standard  rotation  matrix  Rot(^).  As  discussed  in  the  previous  section,  we  assume  that 
v\  is  controllable,  and  V2  is  the  disturbance  input  with  known  bounds.  Safety  is  defined  in  terms 
of  the  unsafe  set: 


Go  =  {qi,  ■  q3}x  {x  <E  X  :  x2r  +  y2  <  52} 

which  encodes  the  fact  that  the  two  aircraft  should  never  come  within  5  nautical  miles  of  each 
other. 

The  solution  of  the  three  mode  conflict  resolution  example  is  presented  in  Figure  6.  Each  of  the 
twelve  graphs  in  Figure  6  is  a  plot  in  the  ( xr,yr )  axes,  with  each  of  the  four  rows  representing  an 
iteration  of  the  algorithm,  and  the  three  columns  corresponding  to  discrete  states  qi,q2,  and  q3. 
The  computation  is  initialized  with  the  set  Go,  which  is  shown  as  the  dark  shaded  disks.  A  fixed 
point  is  reached  after  three  iterations.  The  dark  shaded  region  in  the  fourth  row,  first  column  of 
Figure  6  is  the  most  interesting:  it  is  a  plot  of  the  unsafe  set  of  states  in  q\ .  As  long  as  the  relative 
position  of  aircraft  2  with  respect  to  aircraft  1  is  not  inside  this  region,  then  there  exists  a  control 
policy  such  that  the  conflict  is  resolved.  If  the  relative  position  is  inside  the  unshaded  (white) 
region,  then  switching  instantaneously  will  lead  to  a  conflict;  the  aircraft  must  remain  in  q\  until 
the  relative  state  enters  the  light  shaded  region,  then  either  switch  to  q2 ,  or  remain  in  q\. 


Aircraft  landing  example 

The  autopilots  of  modern  jets  are  highly  automated  systems  which  assist  the  pilot  in  constructing 
and  flying  four-dimensional  trajectories,  as  well  as  altering  these  trajectories  online  in  response  to 
air  traffic  control  directives.  The  autopilot  typically  controls  the  throttle  input  and  the  vertical 
and  lateral  trajectories  of  the  aircraft  to  automatically  perform  such  functions  as:  acquiring  a 
specified  altitude  and  then  leveling,  holding  a  specified  altitude,  acquiring  a  specified  vertical  climb 
or  descend  rate,  automatic  vertical  or  lateral  navigation  between  specified  way  points,  or  holding  a 
specified  throttle  value.  The  combination  of  these  throttle- vertical-lateral  modes  is  referred  to  as  the 
flight  mode  of  the  aircraft.  A  typical  autopilot  has  several  hundred  flight  modes  -  it  is  interesting 
to  note  that  these  flight  modes  were  designed  to  automate  the  way  pilots  fly  aircraft  manually: 
by  controlling  the  lateral  and  vertical  states  of  the  aircraft  to  set  points  for  fixed  periods  of  time, 
pilots  simplify  the  complex  task  of  flying  an  aircraft.  Those  autopilot  functions  which  are  specific  to 
aircraft  landing  are  among  the  most  safety  critical,  as  it  is  extremely  difficult,  if  not  impossible,  for 
the  pilot  to  safely  take  over  control  of  the  aircraft  when  the  aircraft  is  close  to  the  ground.  Thus,  the 
need  for  automation  designs  which  guarantee  safe  operation  of  the  aircraft  has  become  paramount. 
Testing  and  simulation  may  overlook  trajectories  to  unsafe  states:  “automation  surprises”  have 
been  extensively  studied  [18]  after  the  unsafe  situation  occurs,  and  “band-aids”  are  added  to  the 
design  to  ensure  the  same  problem  does  not  occur  again.  We  believe  that  the  ability  to  compute 
accurate  reachable  sets  inside  the  aerodynamic  flight  envelope  under  flight  mode  switching,  may 
help  prevent  the  occurrence  of  automation  surprises. 

A  simple  point  mass  model  for  aircraft  vertical  navigation  is  used,  which  accounts  for  lift  L,  drag 
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D,  thrust  T,  and  gravity  mg  (see  [6]  and  references  therein).  State  variables  are  aircraft  height  z, 
horizontal  position  x,  velocity  V  =  \J x2  +  i2  and  flight  path  angle  7  =  tan_1(|).  Inputs  are  thrust 
T  and  angle  of  attack  a,  where  aircraft  pitch  0  =  7  + a.  The  equations  of  motion  can  be  expressed 
as  follows: 

^  [T  cos  a  —  D(a,V)  —  mg  sin  7] 

•  [T  sin  a  +  L(a,  V)  —  mg  cos  7] 

V  cos  7 

V  sin  7 

The  functions  L(a,  V )  and  D(a,  V )  are  modeled  based  on  empirical  data  [19]  and  Prandtl’s  lifting 
line  theory  [20]: 

L(a,  V )  =  \pSV2CL{a ),  D{a ,  V )  =  \pSV2CD(a ), 

where  p  is  the  density  of  air,  S  is  wing  area,  and  Cl  (a)  and  Cd(o)  are  the  dimensionless  lift  and 
drag  coefficients. 

In  determining  Cl  (a)  we  will  follow  standard  auto-lander  design  and  assume  that  the  aircraft 
switches  between  three  fixed  flap  deflections  6  =  0°,  5  =  25°  and  5  =  50°  (with  slats  either 
extended  or  retracted),  thus  constituting  a  hybrid  system  with  different  nonlinear  dynamics  in 
each  mode.  This  model  is  representative  of  current  aircraft  technology;  for  example,  in  Airbus 
cockpits  the  pilot  uses  a  lever  to  select  among  four  predefined  flap  deflection  settings.  We  assume 
a  linear  form  for  the  lift  coefficient  Cl  (cn)  =  h$  +  4.2a,  where  parameters  ho°  =  0.2,  /1250  =  0.8  and 
/1500  =  1.2  are  determined  from  experimental  data  for  a  DC9-30  [19].  The  value  of  a  at  which  the 
vehicle  stalls  decreases  with  increasing  flap  deflection:  aoiax  =  16°,  =  13°,  al?oax  =  11°;  slat 

deflection  adds  7°  to  the  amax  in  each  mode.  The  left  side  of  Figure  7  gives  a  graphical  summary 
of  the  possible  configurations.  The  drag  coefficient  is  computed  from  the  lift  coefficient  as  [20] 
Co  (a)  =  0.041  +  0.045(7)]  (a)  and  includes  flap  deflection,  slat  extension  and  gear  deployment 
corrections.  So  for  a  DC9-30  landing  at  sea  level  and  for  all  a  €  [—5°,  a™ax],  the  lift  and  drag  terms 
in  (4)  are  given  by 

L(a,  V )  =  68.6  (h5  +  4.2 a)V2  D(a,  V )  =  (2.81  +  3.09  (h5  +  4.2a)2)F2 

Flap  deflection  dynamics  model:  In  reality,  the  decision  to  move  from  one  deflection  setting  to 
another  can  occur  at  any  time,  but  approximately  10  seconds  are  required  for  a  25°  degree  change 
in  flap  deflection.  For  our  preliminary  implementation,  we  have  chosen  to  ignore  the  continuous 
dynamics  associated  with  discrete  mode  switching,  allowing  the  flaps  and  slats  to  move  instantly 
to  their  commanded  positions.  However,  if  such  instantaneous  controlled  switches  were  always 
enabled  then  the  system  would  be  zeno;  therefore,  we  introduce  Delay,  or  transition  modes  Ot,  25f 
and  50t,  which  use  the  envelopes  and  flight  dynamics  of  the  regular  modes  0 u,  25 d  and  50 d  and 
include  a  timer  to  model  the  actual  delay  in  flap  change  (the  discrete  automaton  is  shown  on  the 
right  side  of  Figure  7).  A  regular  mode  may  make  a  controlled  switch  to  a  transition  mode,  so 
flight  dynamics  can  be  changed  instantly.  Transition  modes  have  only  a  timed  switch  at  t  =  tdeiay; 
so  controlled  switches  will  be  separated  by  at  least  fdelay  time  units  and  the  system  is  nonzeno.  For 
the  executions  shown  below,  fdelay  =  0.5  seconds. 

Landing:  The  aircraft  enters  its  final  stage  of  landing  close  to  50  feet  above  ground  level  ([19,  21]). 
Restrictions  on  the  flight  path  angle,  aircraft  velocity  and  touchdown  (TD)  speed  are  used  to 
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determine  the  initial  safe  set  Eg: 


'  z<  0 

v  >  V7ta11 
<  V  <  vmax 

V  sin  7  >  zo 
.  7  <  0 


landing  or  has  landed 
faster  than  stall  speed 
slower  than  limit  speed 
limited  TD  speed 
monotonic  descent 


u  < 


z  >  0 

V  >  V7taU 

V  <  Vmax 


7>  -3° 
,  7  <  0 


aircraft  in  the  air 
faster  than  stall  speed 
slower  than  limit  speed 
limited  descent  flight  path 
monotonic  descent 


(5) 


We  again  draw  on  numerical  values  for  a  DC9-30  [19]:  stall  speeds  Vg^a11  =  78  m/s,  =  61  m/s, 
V/j]/11  =  58  nr/s,  maximal  touchdown  speed  zq  =  0.9144  m/s,  and  maximal  velocity  Hmax  =  83 
m/s.  For  passenger  comfort,  the  aircraft’s  input  range  is  restricted  to  T  G  [0  kN,  160  kN]  and 
a  G  [0°,  10°]. 

The  interior  of  the  surface  shown  in  the  first  row  of  Figure  8  represents  Eo  for  each  of  the  Ow,  25 d 
and  50d  modes.  The  second  row  of  the  figure  shows  the  safe  envelope  E  when  there  is  no  mode 
switching.  Portions  of  Eo  are  excluded  from  E  for  two  reasons.  States  near  z  =  0  correspond  to  low 
altitudes  and  are  too  close  to  the  ground  at  steep  flight  path  angles  to  allow  control  inputs  time  to 
prevent  the  plane  from  crashing.  States  close  to  the  stall  velocity  correspond  to  low  speeds  where 
there  is  insufficient  lift  and  the  flight  path  angle  becomes  steeper  than  that  allowed  by  the  flight 
envelope.  This  latter  condition  holds  throughout  the  very  narrow  range  of  speeds  allowed  in  mode 
On,  with  the  result  that  only  post-touchdown  states  (z  <  0)  are  controllable  in  this  mode.  The 
third  row  shows  how  E  can  be  increased  if  switches  are  permitted  (for  example,  mode  On  becomes 
completely  controllable).  Mode  50ci  is  the  best  to  be  in  for  landing  and  there  is  no  difference  in  E 
with  or  without  switching  enabled.  The  fourth  row  shows  slices  of  the  set  in  the  third  row,  taken 
at  z  =  3  meters.  The  light  grey  regions  are  unsafe  G  and  the  dark  grey  are  safe  E.  The  figure 
shows  that  modes  On  and  25 d  are  safe  only  because  there  exists  a  discrete  switch  to  a  safe  state  in 
another  mode. 


4.7  Overapproximations  of  Reachable  Sets 

In  the  reachable  set  computation  of  the  previous  section,  computational  complexity  introduced 
by  increasing  the  dimension  of  the  continuous  state  is  a  limiting  factor:  the  technique  suffers 
from  exponential  growth  with  respect  to  the  continuous  dimension.  To  challenge  this  curse  of 
dimensionality,  we  are  investigating  techniques  for  computing  overapproximations  to  the  reachable 
set  of  states. 

In  [22],  we  devise  and  implement  a  method  based  on  projecting  the  reachable  set  of  a  high  dimen¬ 
sional  system  into  a  collection  of  lower  dimensional  subspaces  where  computational  costs  are  much 
lower.  We  formulate  a  method  to  evolve  the  lower  dimensional  reachable  sets  such  that  they  are 
each  an  overapproximation  of  the  full  reachable  set,  and  thus  their  intersection  will  also  be  an  over¬ 
approximation  of  the  reachable  set:  the  method  uses  a  set  of  lower  dimensional  Hamilton- Jacobi 
PDEs  for  which,  in  any  projection,  the  set  of  disturbance  inputs  is  augmented  with  the  unmodeled 
dimensions.  For  the  three  dimensional  aircraft  collision  avoidance  problem  of  the  previous  section, 
the  result  of  performing  the  computation  in  the  ( xr,yr )  space,  and  then  back-projecting  this  set 
to  form  a  cylindrical  overapproximation  of  the  reachable  set  in  three  dimensions,  is  computed  in 
about  40  seconds. 

In  addition,  we  have  studied  polytopic  overapproximations.  Here,  the  main  idea  is  to  propagate 
a  set  of  supporting  hyperplanes  of  the  initial  set  such  that,  at  any  time  instant,  the  bounded 
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intersection  of  these  hyperplanes  (known  as  a  polytope)  includes  the  reachable  set.  The  faces  of 
the  polytope  constitute  a  level  set  of  the  special  polytopic  level  set  function.  This  approach  can  be 
successfully  applied  to  feedback  linearizable  nonlinear  systems,  linear  dynamic  games,  and  systems 
that  can  be  represented  as  linear  systems  with  norm  bounded  perturbations  [23].  Reduces  the 
three  dimensional  collision  avoidance  problem  computation  time  from  15  minutes  (time  to  compute 
exact  set)  to  about  1  second  [23].  To  the  best  of  our  knowledge,  the  proposed  method  computes 
an  over-approximation  of  the  unsafe  set  faster  than  any  other  known  method  [24]. 

A  second  technique  is  to  overapproximate  nonlinear  dynamics  by  linear  dynamics  with  bounded 
disturbance,  and  to  overapproximate  reachable  sets  by  polyhedra.  Propagation  of  sets  is  then 
simply  described  by  the  propagation,  through  linear  dynamics,  of  linear  inequalities  describing 
the  polyhedral  set.  With  this  technique  however,  comes  the  problem  of  the  potentially  unlimited 
number  of  inequalities  needed  to  describe  the  polyhedral  estimate.  This  problem  has  been  addressed 
in  the  literature,  yet  previous  solutions  fall  short  in  several  ways:  exact  pruning  of  polytopes  [25] 
requires  complex  data  structures  to  be  maintained  and  the  number  of  “active”  hyperplanes  could 
still  be  unlimited;  the  use  of  simple  polyhedral  shapes  for  bounding  the  polyhedral  set  [26]  is  itself 
computationally  intensive;  the  ellipsoidal  overapproximation  to  the  reachable  set  [27,  28,  29],  while 
an  elegant  and  efficient  representation  of  the  set,  is  generally  too  conservative  an  overapproximation. 

We  have  developed  a  technique  which  combines  the  efficiency  of  the  ellipsoidal  data  structure 
with  the  tightness  of  polyhedral  overapproximation.  The  heart  of  our  technique  is  a  novel  pruning 
procedure  based  on  the  Maximum  Volume  Ellipsoid  (MVE)  contained  in  the  polyhedron.  The  basic 
idea  of  the  method  is  to  rank  the  constraints  by  their  distance  (measured  in  the  norm  induced  by 
the  ellipsoid)  to  the  center  of  the  MVE,  and  delete  the  ones  that  appear  to  be  least  relevant. 
Conceptually,  constraints  that  are  far  away  from  the  center  are  less  relevant  than  those  that  are 
close.  This  approach  will  allow  us  to  perform  guaranteed  pruning  (where  only  provably  redundant 
constraints  are  deleted). 

For  example,  consider  a  discrete-time  linear  system  (which  represents  a  linear  overapproximation 
of  the  nonlinear  systems  of  the  previous  section): 

x(t  +  1)  =  Ax(t)  +  Bw(t ) 

where  x(t)  €  R™  is  the  state  vector,  w(t)  6  Rrf  is  a  bounded  disturbance,  and  A  and  B  are  matrices 
with  consistent  dimensions.  We  assume  that  w(t)  €  W(i)  for  some  bounded  convex  set  W(t). 
Assume  that  x(t)  is  known  to  lie  in  a  polyhedron  V(t),  where  V  C  R”  is  the  intersection  of  a  finite 
number  of  half-spaces, 


v = n  Hi 

i=  1 

Using  the  parameterization  of  half-spaces  introduced  above,  we  can  represent  V  in  the  form 

V  =  {x\  fjx  <  g%,  i  =  1, . . .  ,  m}  =  {x  |  Fx  A  g}  (6) 

where  A  denotes  componentwise  inequality.  Now,  an  ellipsoid  is  defined  as  the  image  of  a  unit  ball 
under  an  affine  transformation 


£  =  {Pu  +  q  |  |H|  <  1} 


(7) 


We  assume  that  A  is  invertible,  which  is  always  the  case  if  the  dynamics  is  discretized  from  some 
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continuous-time  system,  and  we  propagate  the  set  V(t)  through  the  dynamic  matrix  A: 

AP(t)  =  {.x  |  fi(t)T A~lx  <  gi(t),  i  =  1, . . .  ,m(t)} 

=  {x  |  F(t)y4_1x  <  g(t)}. 

To  account  for  the  bounded  disturbance,  we  need  to  find  the  polytope  AV(t)  +  BW(t).  An  easily 
computed  outer  approximation  of  it  is  given  by 

Vt+i  =  {z\fi(t)TA~1z<gi(t)  +  \\fi(t)TA~1B\\i,i  =  li,,.,m(t)} 

=  {z\F(t)z  <g(t)}. 

As  V  evolves  through  hybrid  dynamics,  its  representation  grows  in  dimension.  To  limit  the  com¬ 
putational  requirements  for  storing  and  operating  on  P,  it  is  therefore  necessary  to  occasionally 
simplify  the  description  of  V  by  dropping,  or  discarding,  constraints.  We  will  call  this  procedure 
polyhedral  pruning. 

The  maximum  volume  inscribed  ellipsoid  (MVE) 

We  are  interested  in  finding  the  maximum  volume  ellipsoid  that  is  contained  in  a  polyhedron.  Let 
£  be  an  ellipsoid  of  the  form  (7),  and  V  be  a  polyhedron  parameterized  as  in  (6).  Then,  £  C  V  if 
and  only  if 

SUp  fj X  <  gt,  i  =  1,...  ,  777, 
x€£ 

Using  the  parameterization  of  the  ellipsoid,  this  condition  can  be  expressed  as 

sup  fj  Pu  +  fj  q  <  gu  i  =  1,  •  •  •  ,  m 
||u||2<l 

Maximization  of  the  left  hand  side  expression  gives  the  equivalent  containment  condition 

\\Pfih  +  f?Q  —  9ii  i  =  1,  •  •  •  ,m 


Hence,  we  can  find  the  MVE  inside  a  polyhedron  by  solving  the  convex  optimization  problem 

maximize  log  detP 
subject  to  P  =  PT  y  0 

\\Pfih  +  f?q  <  9i,  i  = 


This  problem  follows  in  the  general  class  of  rnaxdet  problems,  which  can  be  solved  globally  and 
efficiently  using  the  algorithms  described  in  [30].  Moreover,  specialized  algorithms  for  computing 
the  MVE  have  been  proposed  in  [31,  32],  The  MVE  gives  an  optimal  (in  terms  of  volume)  ellip¬ 
soidal  approximation  of  the  polyhedron.  Clearly,  the  MVE  also  provides  an  inner  bound  on  the 
polyhedron,  since  £mve  C  V .  More  surprisingly,  perhaps,  is  that  the  MVE  also  provides  an  outer 
bound  of  the  polyhedron:  the  ellipsoid  obtained  by  scaling  the  MVE  a  factor  n  around  its  center 
is  guaranteed  to  contain  the  polyhedron: 


Anve  U  P  V  gmve  +  n  (Alive  9mve) 


This  is  a  dual  result  to  the  celebrated  Lowner-John  ellipsoid  [33,  34]  (which  treats  the  minimum 
volume  ellipsoid  that  contains  a  convex  set).  Thus,  the  MVE  gives  a  good  approximation  of  a 
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polyhedron,  and  provides  simultaneous  inner  and  outer  bounds. 

Polyhedral  pruning  using  the  MVE 

The  MVE  also  gives  a  natural  measure  for  ranking  the  relevance  of  the  constraints  that  define  V. 
For  a  halfspace  TLi  =  {x  \  fjx  <  g%}.  we  define  the  relevance  ranking  rg  as  the  largest  factor  by 
which  the  ellipsoid  can  be  enlarged  and  yet  still  lie  within  the  halfspace: 

Vi  =  maxja  I  Qmve  V  u(£mve  —  Qmve )  1=  T~ti} 

It  is  easy  to  verify  that  the  relevance  ranking  is  given  by 

Vi  =  idi  -  fIq)/\\Pfih 

This  is,  in  fact,  nothing  more  than  the  minimum  distance  from  the  halfspace  to  the  center  of  the 
ellipsoid  in  the  P 2  metric.  The  relevance  ranking  has  the  following  properties: 

1.  rg  >  1  for  all  i 

2.  Vi  =  1  =>■  fjx  <  9i  active 

3.  Vi  >  n  fjx  <  gi  redundant 

We  define  a  constraint  to  be  redundant  if  the  polyhedron  V  does  not  change  when  the  constraint 
is  removed.  We  will  say  that  the  polyhedron  is  minimal  if  it  does  not  contain  any  redundant 
constraints,  and  we  will  say  that  a  constraint  is  active  if  fj x*  =  gi  for  some  x*  G  V.  A  constraint 
is  (strongly)  redundant  if  and  only  it  is  inactive  for  all  x  G  V.  The  last  property  allows  us  to 
discard  redundant  constraints.  There  is  a  gray  zone  Vi  £  (l,n)  f°r  which  the  relevance  ranking 
neither  proves  that  a  constraint  is  active  nor  that  it  is  redundant.  In  practice,  guaranteed  pruning 
can  often  be  too  conservative  in  that  it  fails  to  discard  many  redundant  constraints.  We  therefore 
use  the  ranking  in  a  heuristic  pruning  procedure:  given  V  =  {x  \  fjx  <  g\  i  =  1 ,  we 

compute  £mve,  sort  the  inequalities  by  their  relevance  ranking,  rg ,  and  keep  the  k  most  relevant 
(the  ones  with  smallest  Vi)-  This  yields  an  approximate  simplified  polyhedron  V,  with  V  C  V.  We 
believe  this  algorithm  to  be  an  effective  method  for  computing  fast  polyhedral  overapproximations 
of  reachable  sets,  in  very  high  dimensions. 

We  have  also  developed  two  methods  for  computing  fast  overapproximations  to  reachable  sets  for 
verifying  safety  properties  of  continuous  and  hybrid  systems.  These  approximative  methods  will 
reduce  computation  times  by  orders  of  magnitude,  and  thus  allow  us  to  treat  problems  of  higher 
dimension  than  before. 

4.8  Summary 

We  have  presented  and  numerically  validated  a  tool  for  determining  accurate  approximations  of 
reachable  sets  for  hybrid  systems  with  nonlinear  continuous  dynamics  and  adversarial  continuous 
and  discrete  inputs.  By  developing  convergent  approximations  of  such  complex  systems,  we  will  be 
better  able  to  synthesize  aggressive  but  safe  controllers.  As  an  example,  the  six  mode  auto-lander 
shows  that  for  envelope  protection  purposes  the  safest  control  decisions  are  to  switch  directly  to 
full  flap  deflection,  but  to  maintain  airspeed  until  touchdown.  With  the  summary  data  from  the 
reachability  analysis,  such  decisions  can  be  made  based  on  local  state  information;  without  it  the 
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auto-lander  may  not  detect  that  low  speeds — while  still  within  the  flight  envelope — lead  inevitably 
to  unsafe  flight  path  angles.  In  this  arena,  we  have  made  real  progress  on  problems  in  a  few 
dimensions,  where  we  have  close  to  exact  solutions.  We  show  that  the  computation  may  be  made 
faster  as  we  loosen  the  numerical  approximation,  however,  even  with  such  approximation,  these 
methods  will  not  scale  beyond  five  or  six  continuous  dimensions. 

In  response  to  this  problem,  we  are  developing  methods  that  gracefully  extend  to  problems  with 
a  much  higher  dimension  of  the  continuous  state  (these  methods  have  been  shown  to  work  well 
in  dimension  up  to  100).  While  these  are  approximate  methods,  they  are,  however,  non- heuristic: 
when  they  determine  that  a  reachable  set  does  not  intersect  some  other  set,  the  result  is  certain. 

Finally,  in  order  to  test  our  reachability  tools,  we  are  designing  and  building  a  system  of  two 
Unmanned  Aerial  Vehicle  (UAVs),  with  on-board  sensing,  navigation,  and  control  systems  enabling 
automatic  flight.  The  Stanford  DragonFly  UAV  Project  currently  consists  of  two  10-foot  wingspan, 
fixed  wing,  unmanned  aircraft,  each  equipped  with  GPS,  inertial  sensors,  and  on-board  computers 
running  QNX  Neutrino.  These  will  be  described  in  an  ensuing  section. 


5  Research  Areas  3  and  4:  Design,  Implementation,  and  Valida¬ 
tion  on  Stanford  DragonFly  testbed  and  on  Boeing  SEC  Plat¬ 
form  (F-15,  T-33) 

In  order  to  validate  the  algorithms  described  above,  we  computed  the  reachable  set  for  a  two- 
aircraft  conflict,  in  which  one  aircraft,  the  “blunderer”,  is  assumed  to  fly  into  the  path  of  the  other 
aircraft,  the  “evader”.  Here,  we  assume  that  the  evader  can  perform  one  of  only  three  “emergency 
escape  maneuvers  (EEMs)”:  this  gives  a  more  conservative  reachable  set  than  if  the  evader  were 
given  full  range  of  maneuvers,  yet  it  provides  a  realistic  and  practical  test  for  collision  avoidance. 
The  reachable  set,  called  the  “blunder  zone”,  is  computed  for  the  three  EEMS:  the  maneuvers 
are  shown  in  Figure  9.  This  example  was  motivated  by  the  problems  of  Closely  Spaced  Parallel 
Approaches  (CSPA)  of  civilian  aircraft  into  airports  with  closely  spaced  runways. 

5.1  Implementation  and  Validation 

The  schematic  for  the  implementation  of  the  new  algorithm  is  shown  in  Figure  10.  The  algorithm 
runs  on  the  ‘own  aircraft’.  It  takes  speed  and  heading  information  from  both  aircraft  and  computes 
the  danger  zone. 

With  the  position  information  from  both  aircraft,  it  then  performs  the  check  of  whether  the  bound¬ 
ary  of  the  danger  zone  encroaches  the  position  of  the  own  aircraft.  If  so,  it  activates  the  necessary 
EEM.  During  the  EEM,  speed  and  heading  information  of  the  blundering  aircraft  will  not  be 
required  unless  they  change  beyond  the  limits  expected. 

Some  issues  and  considerations  regarding  the  implementation  of  the  algorithm  are  as  follows. 

•  Communication  dropouts:  It  is  reasonable  to  expect  communication  dropouts  over  a 
wireless  link.  The  algorithm  takes  this  into  account  and  guarantees  safety  despite  dropouts. 

•  Simplification:  As  it  is,  the  danger  zone  boundary,  for  the  most  part,  is  given  in  the  form 
of  a  series  of  parameterized  explicit  analytical  expressions  (only  one  case  involves  a  bisection 
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search  algorithm).  To  simplify  the  computation  further,  a  slightly  conservative  approximation 
by  a  polygon  is  used,  as  shown  in  Figure  11. 

Using  Microsoft  Visual  C++  ver  6.0  on  a  Dell  Pentium  3  machine,  each  computation  of  this 
approximation  of  the  danger  zone  can  run  in  0.05  sec. 

•  Sensor  inaccuracy:  Sensor  inaccuracy  is  accounted  for  assuming  WAAS-like  characteristics. 
This  is  achieved  by  ‘pushing  out’  the  sides  of  the  polygon  as  shown  in  Figure  11. 

•  Cockpit  display:  A  cockpit  display  of  the  danger  zones  to  the  pilot  have  been  developed 
by  Jennings  [35],  as  shown  in  Figure  12. 

The  danger  zone  computation  is  based  on  a  2-dimensional  kinematic  model  which  is  engineered  to 
contain  realistic  aircraft  trajectories.  These  constitute  approximate  aircraft  models.  A  validation  of 
the  implementation  of  the  algorithm  is  thus  required.  Also,  it  is  useful  in  verifying  the  algorithm’s 
strengths. 

The  implementation  of  the  algorithm  has  been  validated  through  1,106,091  simulations  using  high 
dimensional  dynamic  aircraft  models  representing  a  Piper  Saratoga  as  blunderer  and  a  Cessna 
Caravan  as  evader.  This  choice  of  aircraft  is  because  of  the  collaboration  work  with  Jennings 
[35]  who  performed  piloted  simulation  studies  of  various  displays  (including  these  danger  zones) 
using  these  aircraft,  although  the  motivation  of  the  algorithm  is  towards  commercial  aircraft.  The 
simulation  is  performed  in  3-D  using  non-linear  dynamic  models  whose  basic  control  inputs  are 
desired  bank  angle,  lift  coefficient  and  thrust.  The  algorithm  is  run  on-line  on  the  Cessna  Caravan. 
It  computes  the  danger  zone  boundary  around  the  Piper  Saratoga  and  checks  the  Cessna  Caravan’s 
position  with  respect  to  the  boundary.  Upon  encroachment,  the  Cessna  Caravan  is  made  to  evade 
according  to  the  part  of  the  boundary  that  encroaches  it.  The  horizontal  separation  between  the 
two  aircraft  at  the  closest  point  of  approach,  together  with  the  result  of  each  run,  is  recorded. 

The  objective  of  the  simulation  runs  is  to  validate  the  algorithm  with  realistic  aircraft  models 
under  realistic  conditions  involving  turbulence,  sensor  noise  and  ADS(B)  drop  outs.  While  in 
the  simulation  it  is  necessary  to  choose  only  a  finite  number  of  blunder  trajectories,  recall  that 
the  analysis  considered  all  blunder  trajectories  permitted  by  the  model  and  bounds  on  control 
inputs.  The  validation  approach  is  adapted  from  Honeywell’s  work  [36].  As  presented  in  the 
Honeywell  approach  [36,  37],  three  metrics  are  used  to  evaluate  the  validation  results.  They  are 
the  probabilities  of 

1.  Successful  Alerts  (SA):  this  refers  to  events  where  alerts  are  given,  the  EEMs  are  executed, 
and  no  collisions  occur;  and  if  the  EEMs  were  not  executed,  collisions  would  have  occurred. 

2.  Unnecessary  Alerts  (UA):  this  refers  to  events  where  alerts  are  given,  the  EEMs  are  executed, 
and  no  collisions  occur;  and  if  the  EEMs  were  not  executed,  no  collision  would  have  occurred 
either.  If  the  alert  were  given  because  the  other  aircraft  is  blundering,  the  alert  will  not 
be  deemed  as  unnecessary.  However,  such  cases  are  still  categorized  as  UAs  in  order  to  be 
consistent  with  the  literature  [36,  37]. 

3.  Collision,  Given  an  Alert  (CGA):  this  refers  to  events  where  alerts  are  given,  the  EEMs  are 
executed,  and  yet,  collisions  occur;  and  if  the  EEMs  were  not  executed,  no  collisions  would 
have  occurred. 

The  results,  together  with  those  of  the  Honeywell  AILS  Segmented  [36]  and  MIT  Probabilistic  [38] 
Alerting  Algorithms,  and  the  corresponding  formulae  for  the  metrics  are  given  in  Table  1.  In  the 
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Table  1:  Metrics  for  tl 

iree  algorithms:  %  of  runs. 

750ft 

1700ft 

2500ft 

Formula 

CSPA 

CSPA 

MIT 

CSPA 

AILS 

MIT 

SA 

CD+FA 

CD+FA+MD+IC+LA 

100% 

100% 

99.5% 

100% 

99.9% 

99.3% 

UA 

1C+FA 

CD+FA+MD+IC+LA 

98.78% 

77.14% 

55.5% 

69.02% 

93.7% 

48.7% 

CGA 

IC+LA 

CD+FA+IC+LA 

0% 

0% 

0.5% 

0% 

0.1% 

0.7% 

table,  the  acronym  CSPA  refers  to  the  new  algorithm  demonstrated  in  this  paper,  and  is  referred  to 
as  the  CSPA  algorithm,  while  AILS  refers  to  the  Honeywell  algorithm  and  MIT  refers  to  the  MIT 
algorithm.  To  compare  the  results  with  the  other  two  algorithms,  some  background  is  first  given. 
The  Honeywell  algorithm  is  designed  for  independent  approaches  with  runway  spacings  of  2500  ft 
or  greater.  Thus,  only  the  results  for  a  runway  spacing  of  2500  ft  are  available.  The  thresholds  used 
in  the  Honeywell  algorithm  to  trigger  alerts  are  obtained  empirically  through  simulating  a  large  of 
number  of  parallel  approaches  with  blunders.  The  thresholds  picked  are  the  ones  that  satisfactorily 
prevent  collision.  It  is  not  stated  in  the  Honeywell  paper  [36]  what  ‘satisfactory’  quantitatively 
corresponds  to.  The  MIT  algorithm,  on  the  other  hand,  is  designed  so  that  the  probability  of 
collision  in  the  event  of  a  blunder  is  0.001.  This  criterion  is  considered  to  be  the  acceptable  level 
of  risk  in  the  Precision  Runway  Monitor  System  studies  [39]. 

The  CSPA  algorithm  has  100%  SA  and  0%  CGA  whereas  both  the  Honeywell  and  MIT  algorithms 
have  some  Collisions,  Given  Alerts  (CGAs).  With  regard  to  Unnecessary  Alerts  (UA),  the  CSPA 
algorithm  has  a  lower  UA  rate  than  the  Honeywell  algorithm  but  a  higher  one  than  the  MIT 
algorithm.  This  latter  comparison  may  appear  to  be  a  disadvantage.  However,  as  explained  earlier, 
alerts  in  the  UA  category  would  be  considered  necessary  when  they  are  raised  by  blunders  that 
have  the  potential  to  later  cause  a  collision  should  their  trajectories  be  modified  mid-course.  This 
viewpoint  is  also  expressed  by  Honeywell  in  their  paper  [36].  The  blunders  that  were  correctly 
rejected  (CR)  by  the  CSPA  algorithm  are  those  that  did  not  have  the  potential  as  yet  to  cause  a 
collision.  In  CSPA,  it  is  more  important  to  ensure  that  False  Alarms  (FAs)  do  not  occur  during 
normal  aircraft  approaches  with  normal  flight  technical  error  (FTE).  With  the  CSPA  algorithm,  no 
False  Alarms  will  be  raised  so  long  as  the  aircraft  maintain  their  longitudinal  separation  according 
to  the  MLS. 

The  MIT  algorithm  might  have  produced  better  SA  and  CGA  results  if  it  were  designed  for  a  lower 
collision-probability  value.  However,  it  must  be  noted  that  it  still  relies  on  a  specified  blunder 
trajectory  set,  whereas  the  CSPA  algorithm  does  not. 

In  simulation,  the  new  algorithm  has  been  successfully  demonstrated  to  work  according  to  expecta¬ 
tions  and  is  able  to  guarantee  safety  as  long  as  the  blunder  trajectories  are  captured  by  the  model 
and  are  within  the  bounds  on  the  control  inputs  assumed. 

5.2  Flight  Demonstration 

This  section  describes  the  flight  demonstration  of  the  implementation  of  the  danger  zone  computa¬ 
tion  algorithm  on  the  Stanford  DragonFly  UAVs  and  the  Boeing  testbed  of  F-15  and  T33  aircraft 
and  discusses  the  results. 
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5.2.1  Dual  UAV  Demonstration 


In  this  flight  demonstrations,  DF2  is  designated  the  Own  Ship  which  will  run  the  CSPA  algorithm 
and  conduct  the  evasive  maneuver  when  required.  DF3  is  designated  the  Other  Ship  which  will 
blunder  in  its  ‘final  approach’.  Segment  1  of  the  Danger  Zone  is  not  used.  Segment  2  is  used  in 
its  place.  In  the  demonstration,  both  aircraft  will  conduct  the  closely  spaced  parallel  approach 
autonomously.  DF3  will  be  programmed  to  blunder  in  its  approach  at  a  certain  point.  The 
algorithm  will  run  on  board  DF2  and  will  trigger  an  EEM  according  to  the  danger  zone  computation. 
When  triggered  to  conduct  the  EEM,  DF2  will  fly  the  EEM  autonomously. 

Prior  to  the  flight  demonstration,  a  ground  taxiing  demonstration  is  conducted.  For  the  ground 
taxiing  demonstration,  ground  taxiing  controllers  had  to  be  designed.  In  addition,  a  navigation 
algorithm  which  blends  GPS  measurements  with  IMU  measurements  through  an  Extended  Kalman 
Filter  is  designed.  The  ground  taxiing  phase  is  conducted  prior  to  the  in-flight  phase  so  that  valuable 
experience  in  conducting  the  demonstration  is  gained  without  risking  the  UAVs  in  the  flight  phase. 
The  benefits  of  the  ground  taxiing  phase  are  as  follows: 

•  Testing  Software  Readiness.  The  ground  taxiing  phase  tests  that  the  CSPA  algorithm,  the 
server  and  the  client  processes  can  communicate  with  each  other  and  run  simultaneously  at 
the  desired  rates.  It  also  shows  that  the  CSPA  algorithm  has  been  coded  correctly  to  generate 
the  results  that  match  the  Matlab  version  of  the  algorithm. 

•  Partially  Testing  System  Readiness.  This  phase  also  partially  checks  that  the  dual  UAVs 
are  ready  for  operations.  This  includes  checking  that  the  necessary  logistics  are  available  for 
engine  starting,  battery  charging,  refuelling,  etc,  that  the  UAVs  can  be  switched  between 
manual  and  flight  computer  control,  that  the  installed  sensors  work,  and  that  the  procedures 
for  operations  run  smoothly. 

After  the  ground  taxiing  demonstration,  parameter  identification  of  the  flight  dynamics  of  the 
aircraft  are  obtained  experimentally  and  the  flight  controllers  for  airspeed  hold,  heading  hold,  bank 
angle  hold,  altitude  hold,  and  straight  line  tracking  are  designed  and  flight  tested.  Thereafter, 
the  flight  demonstration  could  commence.  For  the  demonstrations,  each  UAV  has  a  top  level 
architecture  that  is  shown  in  Figure  13.  Below,  the  ground  taxiing  results  are  first  presented 
followed  by  the  flight  results. 


Ground  Taxiing  Demonstration 

Figure  14  shows  the  demonstration  conducted  at  Moffett  Federal  Airfield.  About  10  runs  are  made, 
each  with  a  different  relative  start  position.  In  all  of  the  runs,  the  distance  between  the  two  UAVs 
at  the  closest  point  of  approach  is  greater  than  the  threshold  of  8m.  The  results  of  two  runs  are 
shown  in  Figures  15  and  16.  In  Figure  15,  the  EEM  turns  out  to  be  the  one  associated  with 
Segment  2  of  the  danger  zone,  i.e.,  acceleration  and  turn  to  45  deg.  In  Figure  16,  the  EEM  is  the 
one  associated  with  Segment  3  of  the  danger  zone,  i.e.,  coast  and  turn  to  60  deg.  For  the  former 
case,  the  separation  distance  between  the  two  UAVs  almost  exceeded  the  threshold.  This  is  because 
DF3  turned  at  almost  10  deg/s  in  the  experiment  whereas  a  lower  value  of  5  deg/s  is  assumed  to  be 
its  maximum  turn  rate  in  computing  the  danger  zone.  In  the  latter  case,  the  separation  distance 
between  the  two  UAVs  is  above  the  threshold. 
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The  completion  of  the  ground  demonstration  shows  that  the  CSPA  algorithm  can  run  on  this  dual 
UAV  test  bed  in  real-time. 


In-Flight  Demonstration 

To  reiterate,  the  purpose  of  the  in-flight  demonstration  is  to  show  that  the  algorithm  can  be 
implemented  and  in  particular,  to  demonstrate  that  it  can  run  in  real-time. 

Figure  17  shows  the  demonstration  conducted  over  Moffett  Federal  Airfield.  Over  twenty  runs 
are  made  with  DF3  blundering  and  DF2  performing  the  EEM  as  triggered  by  the  danger  zone 
computation  algorithm.  All  runs  are  successful.  Sample  results  are  shown  in  Figure  18.  The 
former  emulates  an  independent  approach  scenario  where  there  is  no  restriction  in  the  longitudinal 
spacing  between  the  two  aircraft  during  normal  approaches.  The  latter  emulates  a  dependent 
approach  scenario. 

Note  that  the  separation  threshold  is  not  violated.  The  flight  demonstration  shows  the  implemen¬ 
tation  of  the  algorithm  and  in  particular,  demonstrates  that  it  can  run  in  real-time. 


5.3  SEC  Capstone  Demonstration 

The  Stanford  algorithm  of  “dual  vehicle  conflict  resolution”  for  the  Software  Enabled  Control  (SEC) 
Capstone  demonstration  sponsored  by  DARPA  and  Boeing  was  evaluated  successfully  on  the  Boeing 
Open  Control  Platform  (OCP).  The  technology  is  based  on  the  Reachable  Set  Computation  done 
by  Ian  Mitchell  and  Rodney  Teo.  The  experiment  involved  two  aircraft  (F-15  and  T-33),  where  one 
aircraft  is  assigned  as  pursuer  (F-15)  and  the  other  as  evader  (T-33).  The  flight  tests  were  performed 
at  Edward  Airforce  Base  in  June  2004,  and  Stanford  had  six  sorties  during  three  flight  experiments. 
The  algorithm  demonstrated  the  guaranteed  safety  (maintaining  a  minimum  separation  distance) 
for  every  single  flight  test  performed  (there  was  a  single  failure  of  the  algorithm  due  to  an  improper 
initialization  of  the  OCP). 

Figure  19  shows  snapshots  of  a  flight  experiment  displayed  on  the  Experimental  Controller  Display 
installed  on  both  pursuer  and  evader.  As  seen  in  Figure  19(a),  as  pursuer  blunders  into  the  path  of 
the  evader,  the  algorithm  computes  the  corresponding  danger  zone  based  on  position  and  velocity 
information  of  two  aircraft,  yet  the  evader  does  not  commence  an  emergence  escape  maneuver 
(EEM).  As  the  pursuer  enters  a  caution  zone  (which  is  3,000  ft  away  from  the  danger  zone)  as 
shown  in  Figure  19(b),  the  algorithm  alerts  the  pilot  by  displaying  a  warning  message  as  well  as 
changing  the  color  of  the  danger  zone.  Finally,  once  pursuer  touches  the  danger  zone,  the  algorithm 
commences  an  EEM  immediately  based  on  which  part  of  the  danger  zone  the  pursuer  touches. 
Figure  19(d)  shows  the  resulting  separation  distance  between  aircraft  during  this  experiment  and 
it  is  seen  that  the  minimum  separation  distance  (1  mile)  was  well  maintained. 
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Figure  6:  The  rows  indicate  the  iterations  of  the  algorithm,  the  columns  indicate  the  discrete  state. 
Each  figure  displays  a  projection  of  part  of  the  unsafe  set  onto  the  (xr,yr)  axes.  A  fixed  point  is 
reached  after  3  iterations.  Angular  velocity  Wj  =  1  in  q^. 
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Figure  7:  Left:  lift  coefficient  Cx(a)  model  for  the  DC9-30  [19].  Circles  located  at  (a™ax,  C7,(a™ax)) 
indicate  the  stall  angle  and  the  corresponding  lift  coefficient  in  each  mode.  Right:  Discrete  tran¬ 
sition  graph  of  slat  and  flap  settings.  Solid  lines  are  controlled  switches  (<te  in  this  version  of  the 
reachability  problem)  and  dashed  lines  are  uncontrolled  switches  (og). 
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Figure  10:  Schematic  showing  Algorithm  Implementation. 
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Figure  11:  Conservative  Approximation  of  Danger  Zone  for  On-line  Computation. 


Figure  12:  Cockpit  Displays  (for  aircraft  on  the  right)  showing  a  normal  approach  (Courtesy  of 
Chad  W.  Jennings  [35]). 
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Figure  13:  Top  Level  Architecture  of  each  UAV. 


Figure  14:  Ground  Taxiing  Demonstration. 
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Figure  15:  Ground  Taxiing  Demonstration:  Accelerate,  Turn  45  deg  EEM. 
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Figure  16:  Ground  Taxiing  Demonstration:  Coast,  Turn  60  deg  EEM. 
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Figure  17:  In-Flight  Demonstration. 
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Figure  18:  Flight  Demonstration:  Independent  Approach  Scenario. 
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(a)  Pursuer  blunders  into  the  path  of  evader  and 
the  algorithm  computes  corresponding  danger  zone 
displayed  as  blue  line  segments. 


(c)  Pursuer  touches  the  danger  zone  and  the  algo¬ 
rithm  commences  an  EEM. 


(b)  Pursuer  enters  a  caution  zone  and  a  warning 
message  is  displayed  and  the  color  of  the  danger 
zone  turns  into  yellow. 


(d)  Separation  distance  between  aircraft  during  the 
experiment. 


Figure  19:  Snap  shots  of  the  flight  experiment  for  the  Stanford  dual  vehicle  conflict  resolution 
algorithm  on  the  Boeing  SEC  Platform  (T-33,  F-15). 
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